function [Grid,TrMat]=AR1Approx_Rouwenhorst(N,Ybar,Rho,Sig,JumpProb)
% Note: y_t+1 = Y_bar + Rho*y_t +Jump*Eps, where Eps~N(0,Sig)
% Output: 
%   1. Grid:    Grid of State Values
%   2. TrMat:   TrMat_{i,j} denotes the prob. from state j to i.
%% Preliminary
% CHECK IF abs(rho)<1 
if abs(Rho)>=1
    error('The persistence parameter Rho is not smaller than 1!');
end

% CHECK IF n_R IS AN INTEGER GREATER THAN ONE.
if N<1.50001
    error('The number of grid points N must be greater than one.');  
end

% CHECK IF n_R IS AN INTEGER.
if mod(N,1)~=0 
    warning('The number of grid points is not an integer! Rounded value will be used!')
    N=round(N);
end
if nargin<=4
    JumpProb = 1;
end

% CHECK IF n_R IS AN INTEGER.
if mod(N,1)~=0 
    warning('The number of grid points is not an integer! Rounded value will be used!')
    N=round(N);
end
%% Construct the Grid and Transition Matrix
% Step 1: Conditional on Jump
p = (1+Rho)/2;
q = (1+Rho)/2;
psi = sqrt(N-1)*Sig/sqrt(1-Rho^2);
mu = Ybar/(1-Rho);

Grid = mu+psi*linspace(-1,1,N)';
TrMat = [p 1-p;1-q,q];
ii = 2;
while ii<N
    Mat_1 = [TrMat,zeros(ii,1);zeros(1,ii+1)];
    Mat_2 = [zeros(ii,1),TrMat;zeros(1,ii+1)];
    Mat_3 = [zeros(1,ii+1);TrMat,zeros(ii,1)];
    Mat_4 = [zeros(1,ii+1);zeros(ii,1),TrMat];
    TrMat = p*Mat_1+(1-p)*Mat_2+(1-q)*Mat_3+q*Mat_4;
    TrMat(2:end-1,:) = TrMat(2:end-1,:)/2;
    ii = ii+1;
end
TrMat = TrMat';
TrMat = bsxfun(@rdivide,TrMat,sum(TrMat));

% Step 2: Conditional on no Jump
if JumpProb<1
    TrMat = TrMat*JumpProb;
    for ii=1:N
        Gap = bsxfun(@minus,Grid,Ybar+Rho*Grid(ii));
        [~,Ind] = min(abs(Gap),[],1);
        TrMat(Ind,ii) = TrMat(Ind,ii)+(1-JumpProb);
    end
end

